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Abstract 

We revisit the challenging problem of finding an efficient Monte Carlo (MC) 
algorithm solving the constrained evolution equations for the initial-state QCD ra- 
diation. The type of the parton (quark, gluon) and the energy fraction x of the par- 
ton exiting emission chain (entering hard process) are predefined, i.e. constrained 
throughout the evolution. Such a constraint is mandatory for any realistic MC for 
the initial state QCD parton shower. We add one important condition: the MC 
algorithm must not require the a priori knowledge of the full numerical exact solu- 
tions of the evolution equations, as is the case in the popular 'Markovian MC for 
backward evolution" . Our aim is to find at least one solution of this problem that 
would function in practice. Finding such a solution seems to be definitely within 
the reach of the currently available computer CPUs and the sophistication of the 
modern MC techniques. We describe in this work the first example of an efficient 
solution of this kind. Its numerical implementation is still restricted to the pure 
gluon-strahlung. As expected, it is not in the class of the so-called Markovian MCs. 
For this reason we refer to it as belonging to a class of non-Markovian MCs. We 
show that numerical results of our new MC algorithm agree very well (to 0.2%) 
with the results of the other MC program of our own (unconstrained Markovian) 
and another non-MC program QCDnuml6. This provides a proof of the existence of 
the new class of MC techniques, to be exploited in the precision perturbative QCD 
calculations for the Large Hadron Collider. 
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1 Introduction 



The unprecedented experimental precision of the forthcoming experiments at the Large 
Hadron Collider (LHC), in terms of apparatus resolution and event statistics, will have to 
be matched by a far better precision of the theoretical calculations in the strong interac- 
tion sector than available at present. The well established theory of strong interactions, 
Quantum Chromodynamics (QCD), is in principle able to provide very precise predictions 
for the high energy scale (mass, transverse momentum, momentum transfer) processes. 
The perturbative predictions of QCD are obtained within one of two very different cal- 
culational frameworks: the so-called matrix element (ME) calculations and models of the 
parton shower (PS) type implemented in the Monte Carlo (MC) event generators. For a 
more detailed review of these methods, see for example ref . . In the ME calculations 
the basic ingredients are real- and virtual-emission matrix elements evaluated in the fixed- 
order perturbative QCD, for the hard process at the high energy scale, embedded in the 
standard Lorenz- invariant phase space (LIPS). The fixed-order ME is combined with the 
parton distributions (PDFs) describing lower energy multiple emissions in an inclusive 
manner (integrated over the transverse momenta). On the other hand, the PS framework 
offers a fully exclusive picture, down to hadronization energy scale, that is the true MC 
events with explicit 4-momenta, for all multiple soft and collinear emissions associated 
with the hard process - the same emissions as are encapsulated in the PDFs of the ME 
approach. However, the classic PS implements the hard process only at the Born (tree) 
level. 

The above two complementary approaches have their strong and weak points of their 
own. Without entering into details, we may safely say that it is absolutely mandatory 
to combine the virtues of the two approaches if one hopes to ever achieve a significant 
improvement of the precision of the QCD predictions, for a wide class of observables (not 
only total rates); see conclusions of ref. PQ. 

There were numerous attempts to combine ME calculations with the parton shower 
approach beyond the leading order, the most elaborate being the recent one of Frixione 
and Webber [2] . However, none of them are fully satisfactory and there are more proposals 
in this direction; see for instance ref. |3] . There seems to be a growing consensus that part 
of the problem is in the fundamental formulation of the PS models implemented in the 
PS MC. All these models are of the Markovian 1 type, in which the branching process (the 
binary decay of the parton) continues until the boundary of the phase space is hit; the 
number of branchings (emissions) is known at the very end of the branching process. This 
is in stark contrast to the ME approach, where the number of partons involved is defined 
at the very beginning, and the integral over standard LIPS is evaluated for a given ME. 

1 Since the adjective "Markovian" is (ab)used for a wide range of the phenomena, let us state that 
we understand by the Markovian process a walk in a multiparameter space with the consecutive steps 
labelled with the continuous time variable. The rule governing single steps forward ignores the past 
history of the walk. The iterative solution of the QCD evolution equations can be interpreted as a finite 
Markovian process, limited by the maximum time. A Markovian MC implements this process in a natural 
way. In such a MC the number of steps is known at the very end of the MC algorithm. 
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In this sense, the ME approach is basically non-Markovian - this is one of the (principal) 
sources of the difficulties in combining the ME and PS approaches. 

In this paper we do not offer any "silver bullet" solution of the above problems. 
However, we provide one possibly useful cornerstone, in constructing yet another class 
of methods of combining ME and PS methodologies. Our aim is to provide the means 
of reformulating the PS model and the corresponding MC algorithm in a non-Markovian 
way. In fact, we restrict ourselves to an even narrower, but well defined subject of solving 
QCD DGLAP [3] evolution equations using the MC method, which is at the heart of 
any PS MC modelling. We also show, for the initial-state PS (IS PS), that the non- 
Markovian solution of the DGLAP evolution equations emerges in a natural way as an 
alternative solution to yet another long-standing difficulty in the PS MC modelling: the 
problem of the energy constraint. The energy constraint in the IS PS is the requirement 
of constraining to a predefined value the energy of the parton entering the hard process. 
This is so because of a selective nature of the typical hard process ME, typically due to 
narrow resonances. In the typical Markovian PS MC the energy of a parton entering the 
hard process results from many branchings and it is impossible to put any constraint on 
it, in the same way as it is impossible to predefine the number of branchings or the type 
of the parton (quark or gluon) at the end of the branching process. The well known and 
widely adopted work-around is the so-called "Markovian MC for the backward evolution" 
of Sjostrand [5 2 . We shall show that there exists yet another class of MC algorithms with 
the energy constraint, which turns out to be non-Markovian in a natural way. 

Summarizing, the motivation of our search for non-Markovian modelling of the QCD 
evolution equations is that: (a) it is closer to the ME approach, (b) it solves the energy 
constraint problem in the IS PS in a novel way, with potential advantages of its own, as 
discussed below. 

This paper is one of several related works done in parallel, exploiting various aspects 
of the Markovian-type and non-Markovian-type MC solutions of the QCD evolution equa- 
tions. Basic results of the present work were presented in the conference contributions 
quoted in ref. [7j. The earlier work of ref. |S] presents precision MC evaluations of the LL 
QCD evolution equations 3 using an unconstrained Markovian MC. Although the Marko- 
vian calculations of refs. [HI E] are not our main aim, they form a very valuable baseline 
(benchmark) for the constrained non-Markovian calculations, as the ones presented here. 

This work presents the first successful MC algorithm in the constrained non-Markovian 
class, although restricted to the pure gluon-strahlung in the actual numerical implementa- 
tion. Later on, the authors of this paper have found yet another family of MC algorithms, 
in the same important class of non-Markovian constrained MCs, which will be described 
in the forthcoming ref. |10j . and are even more efficient and easier to implement. How- 
ever, at this early stage it makes perfect sense to collect all possible non-Markovian MC 
algorithms for the QCD evolution equations, simply because it is difficult to foresee which 
of them will be most adequate in the future attempts at combining PS and ME calcula- 
tions. In other words, the richer the menu of the different non-Markovian algorithms at 

2 See also ref. UJ. 

3 This work is extended to NLL in ref. [JJ. 
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our disposal, the better. 

The plan of the paper is the following: in the next section we elaborate more on our 
aims and the general framework of our work. In section 3 we formulate in detail several 
examples of the constrained evolution MC algorithms, and present numerical tests of the 
corresponding computer implementations. A short summary concludes the main result. 
The appendix contains the algebra related to the MC method (multibranching) employed 
in section 3. 



2 MC solutions for QCD evolution equations 

As was already said, we are looking for any possibly non-Markovian, MC solution of the 
QCD evolution equations, with the constraint on the final parton type and its x, the 
energy fraction. Needless to say, for a given x, the solutions of the evolution equations 
obtained from the constrained non-Markovian MC will be identical to those obtained from 
unconstrained MC algorithms, or any other non-MC method - the real difference is in 
the efficiency 

The DGLAP evolution equations in QCD, for the quark and gluon distributions in 
the hadron, are derived in QCD using the renormalization group or diagrammatic tech- 
niques jl]. Let us briefly rederive the iterative solution of these equations. We start, as 
usual, from the evolution equations in the standard integro-differential form: 

i 

J x 

where 

f(-)®g(-)(x) = J dxidx 2 5(x — xiX2)f(xi)g(x 2 ) 

and 7kj(t, z) = ^^-P k j(z). Indices i and k = G, q a , % denote gluon, quark and antiquark, 
while the evolution time is t = ln(Q). The differential evolution equation can be turned 
into the integral equation 

t 

e**^D k (t,x)=D k (t ,x) + J dt 1 e*«^YlyU t i>-)®D j (t 1 ,-)(x), 

t j 

where the IR regulator e is introduced: 

? kj (t,z) = -? 6 kk (t,e)S kj 6(l-z)+y%{t,z), (1) 
7%(t,z) = ? kj (t,z)Q(l-z-e) (2) 

and the Sudakov form factor 

Mt,t )= [ dt"? s kk (t',e) 
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appears. The multiple iteration of the above integral equation leads to: 



OO r- 71 p p 

xD k (t, x) = e-* k ^xD k (t , x) + J2 II J dti ^ _ J dZi 



n=l fco...fc n _i i=l 



to 



e -**(t,*n) J dxQ 




Y\_ Z i^%ki^-i Z i) e 



■® ki _ x (ti,ti-i) 



i=l 



x D ko (t , x )5 ( x - x Y[ Z 



1=1 



where k n = k, and the iterative solution is just a series of integrals ready for integra- 
tion/simulation with the MC method. Note that the solution for distributions of parton 
energies xD k (x) is more convenient, because kernels obey the energy sum rules: 



^2 J dz z7i k 



"z) = 0. 



(3) 



More details can be found in refs. jHUHl- 

It is well known ll a that the above iterative solution can be implemented as a Marko- 
vian process with the probability of every single step forward given by the kernel times 
the Sudakov form factor. Formal derivation requires adding the extra integration variable 
t n+ i, t n+ i > t, in every integral; see ref. 0. However, for our present purpose the above 
iterative solution of the evolution equations is the proper starting point. 

In ref. [8 J it was demonstrated that the high-precision Markovian-type MC solution of 
the evolution equations is feasible and it agrees with the non-MC program QCDnuml6 ^2] 
to within 0.2% over a wide range of x and Q. 

Let us still consider one technical point: the choice of the evolution time. The MC 
algorithm will be more efficient if the t-dependence of the strong coupling constant cts{t) 
is absorbed by a suitable redefinition of the evolution time: 



T 



— — - / ah asih), — 
as{tA) Jt A or 



as(t A ) 
a s (t) ' 



(4) 



The choice of t A is arbitrary. For instance, following the one-loop a^\t) = (2ir)/((3o(t — 
hiA )), we may conveniently choose t A such that a^\t A ) = 2ir / ' (3q (e.g. t A — lnA = 1 
and hence t A = ln(eAo)). In such a case r = ln(£ — InAo). The other choice is t A = to, 
where to is the starting point of the evolution. In either case we have 



D k (r,x) = e-' s >^D k (T ,x)+J2 dx Yl U dr t Q(n - n-i) 

n=l " / ° k ...k n -i L i=l ^ T0 ^° 

• n 1 / n \ 

II Zi)e-***-ito'^ D k0 (r , x )5lx - x Q JJ^i) , 

■i=l -I ^ i=l ' 



dzi 



-*fc(r,T n ) 



(5) 
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where k = k n . The kernel 7 and form factor are redefined slightly: 



$fc(r,r ) 



7T 



dr' <?i k (e) = (T-T )r kk (e). 



(6) 
(7) 



In the LL case 7 is completely independent of r^. In the following we shall usually opt 
for tA = ln(eA )) and as(t A )/n = 2//3 . 

3 Constrained non-Markovian MC algorithms 




Figure 1: Graphical representation of the iterative series of eq. 



3.1 Solution types I and II 

What are the general classes of the constrained MC solutions? Let us write once again 
the iterative solution of the evolution equations convoluted with 4 the parton distribution 
Dk {ro,xo) at the low energy scale r and the hard-process matrix element denoted as 
H(x) (see also fig. ^for the illustration): 



cr =^ / dxH k (x)D k (T,x) 
k 

T 1 

//* ^ OO — 7X n n 

dxH k (x) / dx ^ ^2 II / dTi I 

J ° n=0 k n -i...kiko M=l T £ a o 



(8) 



x g-(T-r„)iJ fe 



n 1 / n \ 

U 7 l^M) e- {Tm)Rk ^ slx-xoJlzA D k0 (r ,x ), 

1=1 -I ^ 8=1 ' 



4 For simplicity we include here an iterative solution of the evolution equations for the single initial- 
state hadron, but our real interest is the case with two initial-state hadrons. 
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where we define k = k n and J~[ = 1 in order to keep the formula compact. In the LL 



k=l 



kernels and form factors simplify and from eq. (J3J) it follows that 



5. 



$fc(T,T ) = (r- r )i2 fc 



l-e 







dz z9%(z) 



(9) 



For the purpose of future discussion we define here additional virtual form factors: 



R 







dz z?%{z) , R' k = Y, R jk = Rk - R 



kk- 



(10) 



Let us discuss basic limitations and possible solutions for the MC implementation of 
the above series of multidimensional integrals. As already stated, in the ISR case, since 
there are narrow resonances in the hard-process function H(x), the variable x has to be 
the first one generated in the MC algorithm, i.e. it has to be the outermost integration 
variable. Similarly it is better to keep k = k n as the outermost summation variable as 
well. 

The central issue is the following: How do we treat the variable xq7 There are two 
possible options. In the first option (I) xq is kept as a second outermost integration 
variable, next to x, i.e. it is generated in the MC as a second variable. In the second 
option (II) Xq is treated as one of the last variables in the MC - in fact it is derived from 
the other ones using energy constraint. 

The following formula describes the first case: 



o = f — [ dx H k (x) D kko (r, — tq) D ko (T ,x ) 

ik1 J Xo J v Xo J 

T 1 

OO r- n p p 



n=0 fe n _i...fei L i=l 



n-i 



n -i / n \ 

U^M) e-^-^^-i slz-l[zA 

i=l -I ^ i=l ' 



with k = k n as usual. We shall refer to this option as to a solution type I; this is the 
scenario that was implemented for the QED ISR pure bremsstrahlung in several YFS-type 
MC programs, starting from the prototype of ref. The main technical difficulty is 

the implementation/elimination of the 8(z — Yl z i) function. This problem was solved in 
QED by eliminating the delta function with integration over the z of the hardest photon 6 
(the largest 1 — z). This scenario looks definitely feasible, and the first working example 
will be described in a separate work; see ref. [TU] . 

5 In the NLL case an additional r dependence through as(r) will invalidate such a simple relation. 
6 In this QED case the integration over xq is rather trivial because one starts from D(tq, z) = 5(1 — z). 
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The second scenario, referred to as solution type //relies on the fact that the starting 
parton distribution for typical hadron beam particle D ko (r ,x) can be relatively well (by 
MC standards) approximated by a power-like function D ko (r ,x) ~ x~ l+T < over a wide 
range of ln(x) (i.e. 10~ 4 < x < 1), with the parameter 77 not far from zero. In fact the 
gluon and quark singlet parton distributions of the nucleon at low Q feature rj ~ —0.2. In 
solutions of type II the essential idea is that, in eq. (jHJ), the 5[x — xq YYi=i z i) ls eliminated 
using the integration over xq. It means that xo, contrary to type I, is generated in the MC 
as a last variable instead of a second one. More precisely it is not generated at all, but 
determined as a function of all previously generated variables Xq = x(Yl z i) ; by m eans 
of solving the energy constraint. 

Let us isolate explicitly the small- z limit from the starting parton distribution: 



D ko (ro,x ) = W^(xq) A ko x 



*D 
/>'() 



_ D ko (TQ,X ) 

' ° } = A, x~ l+r > ' 



(12) 



where W^ D is the MC weight to be neglected now and restored later on. Elimination of 
the ^-functions with the help of the xo-integration leads to 



/OO ~ Th n n 

dx H k (x) A ko x~ 1+r > E II / dTi / 

n=0 k n -\...kiko *~ i=l ^ n 

>i — l u 



dzi 



x e -(T-r n )flfc 



n e-^*-^ e ( n Zj - x ) w* o d (x ) 
1=1 j j=i 



(13) 



where x = xf Y\- Zj. For the distributions of quarks and gluons in the proton at a low 
energy scale we have the same i] ~ —0.2. Hence, we may also include all factors z~ v in the 
MC weight W k < 1, to be neglected and restored at a later stage of the MC algorithm 
together with the other details of the D k (r ,xo) function, or in a more sophisticated 
MC algorithm, we may actually generate exactly the distributions z^Tf.^^Zi). In the 
following let us assume the former simpler case with z~ v in the MC weight: 



a 



^2 J dx H k (x) A ko x 



-l+r? 



x e" 



-(r-T„)R k 



G 

k, k. 



E E 

n=0 k n —\...k\ko 



n 



dTi 



dzi 



n 



(14) 



Dk (ro,x ) 
A, t~ 1+7] 



n 



Dk (To,X ) 

A, t~ 1+v 



X 

x 



< h 



where D ko (T ,x ) — > A ko XQ 1+v in the limit x — > 0. 

Let us stress that eq. ()13|) implements the exact iterative distribution of the evolution 
equations. In the following sections we will show the results from the prototype MC 
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based directly on the above expression for the pure bremsstrahlung case. As we shall 
see, it works well for the case of the emission from a quark; however, it has rather low 
acceptance (~ 10~ 3 ) for the emission from a gluon. We shall call such a MC solution, 
directly based on eq. (JT3J), solution type II. a. 

In another method, II. b, we shall reorganize the integration variables in a hierarchical 
way, and use multibranching to isolate the 1/z part of the gluon-to-gluon kernel. Such 
a solution is rather complicated and non-trivial to implement in a general case with 
multiple flavour- changing (quark-gluon) transitions. However, successful implementation 
of the pure gluon-strahlung case, presented in the next section, allows us to claim that 
the efficiency of the MC type II. b is satisfactory, and the gate to practical applications of 
MC solutions of this type is wide open. 



3.2 Constrained Monte Carlo: solutions class II 

As already noticed in ref. jH], in the long emission chain (on average ~ 20 emissions), from 
Q = 1 GeV to Q — 1 TeV, most of the emissions are of bremsstrahlung type, i.e. they 
preserve the identity of the parton on the main line of the chain. It was shown there that 
on the average only about one out of twenty emissions involves flavour transmutation, 
q — > G or G — > q\ the other ones are gluon emissions. With this in mind, it is therefore 
natural to reorganize the iterative solution of the evolution equations in such a way that 
all pure bremsstrahlung adjacent vertices in the emission chain are lumped together into 
segments described by the following universal evolution function: 



d' k (r,z\r ) =e^ T - T ^ Rkk 



T 1 

oo n « « / n 

n=0 i=l T j_ x { V j=1 



(15) 



where the type of parton on the main line, k = Q or k = G, is unchanged. Note that we 
have retained in eq. (fTB^) only part of the virtual form factor R k , namely the Rkk function, 
which matches exactly the real emission kernel of pure bremsstrahlung. The leftover 
R' k is included explicitly in the following eq. (fTHj) . 

The full iterative solution of the previous section can be expressed in terms of the 
product of the above functions and the kernels representing flavour- changing transitions 
q — > G or G — > q in the following way: 



Du(t,x) 



n=0 



k n _ 1 ...,k 1 ,k \i = l 



e /<** e n *> nil ** nH 



kjjtkj_i,j=l,...,n 

xe -(r-r n )K d ' k {r, z '\ Tn ) 



Ti-1 



.1=1 



J=l 



n 



i=i 



+11 

(n~Ti-i)R' k . 



1 d'^n, ^|r<_i) 



(16) 



(n n+1 
x ~ x ° n Zi n 



i=l i=l 



where we employ the usual conventions: k n = k and Yl°j=i = !• 

We say that the above formula implements hierarchical organization of the emission 
chain, because it represents the Markovian process in which each pure bremsstrahlung step 
(segment) in the Markovian random walk (emission chain) is an independent Markovian 
process of its own! This elegant and powerful reorganization of the emission chain is 
proved formally in the separate work of ref . pi] . 

The complete hierarchical formula for the integrated cross section, after elimination 
of the delta function with xq integration, reads as follows: 

r 1 1 

/co n „ n „ ri+l „ 

dx H k ( x ) A ko x~ i+r > j2 J2 n / d ^ n / dz * n / ^ 

n=0 fc n -i-.fci,fc »=1 J »=1 n i=1 n 

fc i #fc J -_l,J=l,.-.,n ri "! u u 



n 



X 



X 



i=l 



(n n+l \ 
i=l i=l ' 



d'ki T h 4\ T i~i) = e 



-(n-n-i)Rk 



En/ ^ / ^ ^(^j?) 5 N - n z $ 

i(i)=0 m=l (4) Q ^ m=l 



(17) 



where 



-(*) 



•7-W = T . -j- 

' — 'u 'o 



(0 



r i-l ; ^0 



X 



< 1. 



(18) 



nLi*nSN" 

The above formula is the starting point in the next two subsections for construction of 
constrained MC algorithms of type II. 



3.2.1 Solution Il.a 

In the straightforward solution of type II, which we call Il.a, the single B-function for both 
flavour- changing emissions and pure bremsstrahlung segments is replaced by the product 
of the individual 6-functions - thus decoupling completely the ^-integration space and 
opening the way for the analytical integrations of the approximate spectra for the purpose 
of the MC generation. More precisely, let us first notice that we are really dealing with 
the single 0-function involving all z variables due to trivial identity 



n+l 

n 

i=i 



\ / n« * s , n« x "1 , n n+l 

i ^m=l n ' ^ m=l ' ^ i=l i=l 



— X 



n+l / n 



i=\ x m=l 



n n *® e n*nn 



n+l n« 



(19) 



r (0 _ 



i=l i=l m=l 
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We are now ready to describe the essence of the MC algorithm of type II. a. We use 
the following identity 

n n+1 raW \ n n+1 j»W 



e(n*nn 

^ i=l i=l m=l 

where the function 



z® - x 



W 



e 



Il.a 



i=i 



n n+1 n 



(0 





7/.a 



e(n*nn 

i=l i=l m=l 



i=l m=l 



X < 1 



(20) 



(21) 



is the Monte Carlo weight. This MC weight will be neglected later on, so that variables 
can be generated according to simplified distributions, and finally the generated events 
will be weighted according to this weight. 

Let us rewrite our master integral (|17|) without any approximations 

r 1 



k X 



a = y^ dx H k (x) A 
k 

x e- (T - r ' l)i? '<(r,x|T n ) 



-1+7] 



e e n 



n=0 k n _ 1 ...,k 1 ,k 

kj^kj_i,j=l,...,n 



1=1 



T»_l 



(ifc(rj,x|rj_i) = e" 



i=i 

nW=0 m=1 (i) m=1 x 



W^xoWZ 



© 



(22) 



where xo = x/( niLi ^ llTi 1 11™= 1 ^ J anc ^ one snou ld remember that (/"-functions pro- 
vide integration over all the zf 1 variables that are implicitly present in the function W® a . 
The enormous advantage of the above procedure is that in the approximate integral we 
can sum up and integrate immediately over all bremsstrahlung segments of the emission 
chain. To see it let us drop the two MC weights W k D Q and W® ' a . Now we can immediately 
sum up and integrate analytically the pure bremsstrahlung subintegrals: 

l 



d'l(ri, x\n-i) = expl (n - t, 



i-i)(-Rkk + I dz ?® k (z) 



(23) 



and hence 



a 



dx H k (x) A ko x- 1+ v J2 II dTi dz 

n=0 fe n _i...,fci,fc i=l J 

n 

ny 8 (z ) g^^-'J^-iU^^!^!^)) 



(24) 



{T-T n )(-R k +Sldz7f k (z)) 



1=1 
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The above looks rather promising, because we are left with the relatively simple problem 
of generating several variables n and Z{ (i < 5 seems sufficient) for the flavour- changing 
emissions, for which the above integrals provide explicit analytical distribution. The MC 
events are attributed with the MC weight Wn. a = W^(x ) Wff a < 1. The key question 
is: What is the acceptance rate for this weight? We did an introductory exercise, imple- 
menting the pure bremsstrahlung version of it. We have found that, unfortunately, the 
acceptance rate for the emission from the gluon line is only about 10 -4 . This inefficiency 
can be traced back to the presence of the 1/z singularity in the G — > G kernel. Namely, 
if we allow for the range x > x m ; n = 10~ 5 we also allow for z to be generated to the same 
low limit. In the case of Tgg^) containing a 1/z part, this creates many events with 
low Z(. Consequently, xo = x/Y\ z i goes very often beyond 1 and the corresponding MC 
weight W® a gets zero value. 

On the other hand, for the quark line the acceptance rate is close to 1, which is clearly 
related to the absence of the component 1/z in < J > qq (z). 

The above numerical exercise indicates that the 1/z part in the ^Pgg(z) has to be 
treated better, as is done in the present case II. a. A more sophisticated treatment of the 
1/z component of the kernel is applied in the solution II. b described in the next section. 

The present solution II. a is still a workable solution. In spite of its very low efficiency, 
due to its relative simplicity, it can still be quite useful for testing other more sophisticated 
solutions. We therefore implemented it also in the MC program, for the moment only in 
the pure bremsstrahlung version. 

3.2.2 Solution Il.b 

In this section we present a solution more sophisticated than the II. a one of the previous 
section: we split the bremsstrahlung kernel 7gg into t wo parts: (B) ~ 1/z and (A) 
= the rest. Then we apply the multibranching 7 , as described in Appendix A, to every 
pure bremsstrahlung segment of the emission chain. Finally, we also treat the O-function 
more selectively than in II. a. The product of the individual 0-functions appears for the 
flavour- changing emissions and part (A) of the bremsstrahlung kernel, while the single 
0-function is left for each segment describing part (B) of the pure bremsstrahlung. Such 
partial decoupling in the ^-integration space still allows for the analytical integration of 
the approximate spectra for the purpose of the early stage of the MC generation. This 
is possible because segments of type (B) in the pure bremsstrahlung parts are integrable 
(to a Bessel function), as shown below, while for the rest we get exponentials in a way 
similar to those in II. a. 

Again, the starting point is the complete hierarchical formula (|17|) for the integrated 
cross section. In the case of the gluonic subintegral d' k=G (Ti, z'^r^i) we reorganize this 
integral to isolate the 1/z part from the kernel. In order to split the , 3 > qq{z) i n two positive 

'Multibranching or multichannelling is the standard MC technique in which the distribution is split 
into a sum of positively defined subdistributions. First the index numbering the distributions is generated. 
Once the subdistribution is chosen, a MC point is generated according to this subdistribution, instead of 
the total distribution; see ref. ^3 f° r details. 
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parts, one of them being 1/z, we have to simplify it first: 

^ G (z) = ?% G {z) w GG (z), PS g (z) = PS G (z) w gg (z), 



e , , _ « n I 1 M ... ,„s _ n ^ ^2 . , ( 25 ) 



Pg G (z) = 2C A ^~ + —J , w GG (z) = (1 - z(l - z)Y < 1. 

Having done that and using the multibranching identity of eq. (j91|) in the appendix we 
may rewrite the gluonic bremsstrahlung subintegral for k = G as follows: 

l 



(26) 

where 

dj*(r i ,Z«|r i _ 1 ) = e-^- 



>< e n / ^ - A) n / d ^ ^(4°) *(* w - n < 

n'( i )=0 fi=l m=l m=l 



(0 



df(r i ,ZW|r i _ 1 ) = e-^'- « 

• (i) n „»(i) „«(i) 



x 



En/ < w « w - ^-i) n / d ^ 5[z^ - n <?> 

n , ( i )=0 W=l m=l ^ m=l 



(27) 



and 



-^ggW — -i 5 Pgg^ z 



1-z 7 ^ v 7 z (28) 
^Sg( z ) = ^ggW' ^ggW = ^gg( z ) ~~ ^gg( z )- 

In the above we used the notation z'^ = {z ,( (\ . . . , z /( ^ {i) ) and z*W = (z*^ , . . . , z,^. 
The MC weight due to the kernel simplification 

W Ga (z®{f®,*P ) )) = n uto^^V")), (29) 

m=l 

depends on the variables z$ after relabelling. A special kind of permutation 2$ — > 
z* , z'f , which we refer to as relabelling, is an important part of the MC algorithm - it 
is defined precisely in the appendix. Since relabelling is just a permutation of z's, we may 
calculate the weight W GG with the z'^\ z'f^ variables before the relabelling: 

n .(i) „/(*) 
W GG =l[w GG (z-®) l[w GG (z'®). (30) 

m=l m=l 
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Until now we made no approximation in our master integral - we only reorganized 
integration variables, in particular isolating the 1/z component in the pure bremsstrahlung 
subintegrals for the gluon emitters. In full analogy to case II. a, in the last step in this 
reorganization we eliminate all variables z[ and a class of 5-functions with the help of the 
identity: 

n+l \, n+l , , n' w s , n n+l 

n/^e(zw-^) n^-n^ © (rhiN-* 

i=l £ i=l ^ m=l ' ^ i=l i=l 

n+l / n n+l ✓ n' (4) \ 

11-11 III II ; ) 

i=l \ i=l i=l V m=l 7 / 

Consequently, from now on we substitute with z[ 

m=l 

On the other hand, we keep 5-functions inside the cIq functions, which will also be treated 
analytically, but separately; see below. 

Now comes the essential step in the algorithm II. b - we define the following MC weight: 




(32) 



n n+l n s n+l n' (i) 

eHiMi^Hnn^ 

j=l j=l * ^ i=l m=l 



X 



(n n+l \ n+l n'^ / s n n+l -> 

n * n zU) - x ) n n @ ( { n n z0) k - - ) 

j=l j=l / i=l m=l ^ ^ j=l j=l J ' 

n+l n' (i) 

= < b (l - ^(Z, Z)) J] II 6 (^m - Z)) , 



i=l m=l 



(33) 



where 

^(^.^....^zo.zw,...,^*")) ^ n.+wor (34 > 

1 lj=i i 1 lj=i 

The new MC weight 

(n n+l , n' (i) , \ 

i=i i=i V m =i / / 



will be neglected later on and restored at the end as a standard MC compensating weight 
8 There are a few other slightly different possible choices of z^ in , which are not discussed here. 
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All these preparatory steps lead us to the following master equation for method II. b, 
still without any approximation, but with clearly defined MC weights and the distributions 
to be generated at the early stage of the MC algorithm: 

T i L1 1 

n+l r dzi i) 



/LXJ tip II n flT-L /» 

dx H k ( x ) A ko x- i+r > J2 E U dr U dz U 



71=0 k n _ 1 ...,k 1 ,k i=l i=l i=l 



i=i - 

x^oXeii-^), 

/(») TJ 



oo n v ' » 

ctf , (r*,A|r i - 1 ) = e-^^ Ell/ 

n '(0=0 m=1 



r (<) 

x n / ^ ^rf) - ^ ^c(z- (i) , z ' (i) ] 



m=l 







(36) 



The functions d^(r,i, Z^\r,i_i) and the variables are really present only for gluon, 
k = G. However, in order to keep the notation compact, we understand that 

^ g (t,,Z(^_0^(1-Z«). (37) 

We also assume Y\n=i = 1> as usua l- The reader should also keep in mind that at this 
stage the integrand of the part d A still depends on the integration variables of d^. 
In the MC algorithm of type II. b all three MC weights are neglected: 

W k D (x ) < b , W GG (z^, z' (i) ) = W u . b (38) 

at the early stage of the MC algorithm and later on generated MC events are weighted 
with Wn.h- Since Wu.b < 1, we may easily transform weighted MC events into unweighted 
events by rejecting some of the MC events in the usual way. 

The MC weight Wn.b was chosen in such a way that once it is neglected, we can 
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perform a lot of analytical integrations: 

T 

/OO Th /l lb /l / L j _1_ /i 

^ ^(x) ^ ^ n / dT * n / n / 



n ^ n+1 



n=0 k n _ 1 ...,k 1 ,k i=l 
kj^kj_i,j=l,...,n 

1 



Ti-l 



1=1 



i=l 



Xe -(r-r n )K df(r,z^ m \r n ) -^rff^Z^IrO 

n r 1 

x]I 4,-^)^ (T, " T '" l)J?tl <C (r*, Air*-!) ^K^VO 



1=1 L 






n+1 


xe(n 








i) = e"^ 





« En/ ^ n / d ^ 

n '(<)=0 m=l m=l g 



= exp I (r, - ( + I dzV® k A (z 



M _ 

i r Q — T i-li 



rfg(r„Z('Vi-l) = e-("-"-'l R "4«(Z» - 1) + 



n*W=l 



(Tj ~ Ti_i)' 



.(») n'W I 



m=l 



•(<) 



m=l 



(39) 



The distribution for the trouble- making component of the kernel Pqq(z) = Pqq(z) 
2Ca/z and Rq G = (a s(t a) /ti)2C 'a, can be calculated analytically: 



^(r^VO =e -(— *-Oflg G 



- 1) 



-.00/ n w n„i /• n \ 

+ w, e ^r^r n / <«* - n *) 

n=l i=l - 70 V 7=1 7 



(40) 



The integral 



n „i / n \ n „ ✓ n 

JJ / dzj <S( Z w - Y[zA = Yl / din* S( lnZ (i) - E ln 

i=i^° V j=1 J i=1 J V 7=1 



^(l/zW)]"- 1 



(41) 
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is just the volume of the simplex and when inserting it in eq. (|4"U|) we find 



<ig(T i ,Z»|r,_ 1 )=e" (r '" T -' |B "° 



ln r ' 



Z (i) ^ n \(n-l)\ v Z« 

71=1 

The above can easily be expressed in terms of the I\ Bessel function: 



(42) 



o*i(2; u) = = 4f J i( 2 v^)- ( 4 3) 

z — ' nlln +1)1 \/u 

n=0 \ / v 

We shall, however, introduce our own notation: 

B'(ri,Z®) = e-*[6(l-Z®)+rioF 1 {2-,-r l ]nZ®)], (44) 
which leads to the simple expression 

d*(T l ,zV\T l _ 1 ) = ^B>((T l -T l _ 1 )RZ G ,Z^, (45) 

which can easily be plugged into a MC program. 

In the above results of the analytical integrations, we easily identify the compact 
analytical expression for the distributions of the z and r variables of the flavour- changing 
emissions (upper layer in the hierarchy). For each pure gluonic segment, there is one 
additional variable Z. 

Since the average multiplicity of the flavour- changing emissions is ~ 1, we may sim- 
ply plug in the integrations over Ti,Zi,Z^ % \ i = l,...,n max into any general-purpose MC 
simulator, for instance into the FOAM program [TBI IT7] . The value of n max = 5 is probably 
more than sufficient for a precision of 10~ 4 and it is feasible for FOAM (up to about 20- 
dimensional distributions), especially because the integrand does not involve any strong 
singularities. Also, generating points according to the higher- dimensional distributions 
will be done very rarely. 

This completes the theoretical description of the constrained MC algorithm of type 
Il.b. 



3.3 Construction of non-Markovian constrained MCs, type II 

In this section we present an actual implementation of some of the constrained MC al- 
gorithms of class II described in the previous sections. Some numerical results are also 
given. 

We shall proceed from simple examples of the MC algorithms for simplified distribu- 
tions, gradually going to more elaborate examples in which the previous, simpler, MC 
examples are used as benchmarks in the numerical tests. It is worthwhile to describe the 
above step-by-step method of creating more and more sophisticated versions of the MC 
algorithm and its numerical realization, because it is an essential part of constructing any 
precision MC event generator, albeit it is rarely explicitly exposed in the literature. It 
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can be of vital interest for any reader interested in the practical aspects of constructing 
MC event generators 9 . 



3.3.1 Benchmark MC for Pqg = 2Ca/z, Poisson-type and inefficient 
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Figure 2: Pure gluon case. Pqg = 2Ca/z. The distribution of x = Yl i Z{ and its ratio to the 
analytical prediction B'(^f,x). 

As a warming-up exercise, let us now work out in detail a MC algorithm calculating 
the following integral, cf. eq. (pfUj): 



J(7) 



dx c?g(r, x\t ) 



f 1 1 

/ dx— B'{j, x), 



(46) 



where 7 = (as^y^C^T — r ) and ej <C 1 
J B'( 7 ,x) = e" 7 



00 „ n „i 



n=l i=l 



(47) 



with the aim of preparing basic tools and setting baseline normalization for the MC 
algorithm of type II. b (similarly as it was done in ref. JB]). On the one hand, the Bessel- 
class function B' is known analytically in terms of a series (JUJ). On the other hand, the 
integral 7( 7 ) can be rewritten as 



— n ,n _■_ 



n=l 



rl dzi 



=1 J £i 



e( l[z j -e 1 

3=1 



(48 



9 Such simplified MC programs existed for many precision MCs for the QED calculations with YFS 
exclusive exponentiation; see for instance ref. 
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remembering that x = rL 2 *- The above integral is easily implementable in the MC, 
which treats the function as a MC weight: W = O^IIJ=i z j ~ e i)- The variable n is 
generated according to the distribution 

J = e~\ I n>0 = ^-L In" -, I = V I n = e-~<e-~^\ (49) 

The variables Zi G (ei, 1) are generated according to the distribution 1/Zi. Once we 
generate suitably long series of MC events (n; Zi, z%, z n ) we calculate the integral using 
the average weight, with the usual expression / = (W)I. In the same MC run we can 
also obtain the distribution B'( r y,x)/x, just by examining the histogram of x = YliZi. 

In the LHS plot of fig. |2 we show the (properly normalized) distribution of x from 
the MC. The acceptance rate ~ 3 x 10~ 4 is rather low - it demonstrates the problem 
with the \ jz component in any MC (also Markovian) in which the starting point of the 
generation of the emission probability is of the Poisson type. This phenomenon is quite 
general. Our numerical example shows the evolution from Q — 1 GeV to Q — 1 TeV. The 
average emission multiplicity in the MC run is about 3.4 for e\ = 10~ 3 . Since the resulting 
distribution of x is known analytically, we can also examine its ratio to the MC result. In 
the RHS plot of fig. |2]we show this ratio (for 1.4 x 10 9 events). It is equal to 1, to within 
the statistical error of order ~ 1%. In the LHS plot we clearly see that the contribution 
~ 5(x) is reproduced by this MC algorithm/program (absent in the analytical program). 

For the purpose of the next exercises we are interested not only in the value of the 
integral, but also in the exclusive distributions. In fig. El we examine the distributions of 
the first four variables = YY n =i an d T i i n the emission chain. The ordered Tj variables 
are generated within the range (r , r) corresponding to Q = 1 GeV and Q = 1 TeV. 
In the following we shall check that the above semi-exclusive distributions are correctly 
reproduced by more sophisticated MC algorithms. In this figure we also include the 
distributions of the emission multiplicity and the MC weight. In the weight distribution 
we exclude zero-weight events. 



3.3.2 Weight-1 algorithm for Pqg = 2CU/ 'z, Bessel type 

The inefficiency of the algorithm described in the previous subsection is mainly due to 
the fact that the emission probability distribution in the integral under consideration is of 
the type P n ~ x n /n\(n — 1)!, Bessel- type for short, while in the MC we actually generate 
a Poisson distribution P n ~ x n /n\ and turn it into a Bessel-type one by the inefficient 
brute-force rejection method. Now we proceed to the next step - we construct a prototype 
algorithm in which a Bessel-type emission probability is used from the start and there is 
no need for the rejection at all. The previous inefficient Poisson-type MC will be useful, 
however, as a precision cross-check for the new one, especially for testing semi-exclusive 
distributions. 

Let us consider almost the same integral 

f 1 1 

J / ( 7 )= / dx-h{x) 1 Fi(2;-7ln(a;)), (50) 
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Figure 3: Pure gluon case. Kernel= 2Ca/z. Distribution of the variables x, = 1X1=1^ anc ' 
of the ordered r< G (r , r). 



which in the multi-integral form looks as follows: 



I'd) 



dx— h(x) 

x 



n „i 



n- — Jn 

n=l i=l J{i 



i=i 



(51) 



where we have removed the unimportant 5(1 — z) component and inserted the test func- 
tion 10 h(x)/x. This integral can be rewritten as 



/V,) / ^h(x) 

X 



7 n ln n - i (l/x) / (n-1)! 



^ //!(//-!)! Vln"- L fl/.r) 11 

n=l v ' N \ 1 / l= \ 



—, — — TT / dlnzi tfflnx - T^ln 



where the internal part of the integrand is conveniently normalized as 

n r i 

I, — LJl 



( 1)\ n pi n 



1 



(52) 
(53) 
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In the following numerical exercises we set it to the constant value h(x) = 1.908359. 
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Figure 4: Pure gluon case. Kernel= iC^jz. Distribution of the variables Xj = [j^ =1 ^ and 
of the ordered T{ G (r , r). 



and we may simulate z variables very easily. Changing variables to i/i = In Zj, we see that 
the distribution in i/i is a uniform distribution over the n — 1 dimensional simplex. There 
are several convenient methods of generating points uniformly within such a simplex. The 
simplest method is to throw randomly n— 1 uniform points Ui G (lnx, 0), % = 1, 2, ... , n — 1 
and order them using any standard method: In a; = u n < u n ^i < ■ ■ ■ < Ui < u = 0. 
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Then we take the differences yi = Ui — i — 1, . . . , n, which by construction fulfil the 
constraint In a; = YTj=iUj- 

The MC algorithm consists of the following steps: first the x is generated according 
to the distribution 

p(x) = - h(x) 7 Q Fi(2; -7 ln(a;)). (54) 
x 

Next the number of emissions n is generated according to a (normalized) Bessel-type 
probability 11 distribution 12 : 

/ \- 1 'y n \n n ~~ 1 (l/x) ^ 

P„=(7oF 1 (2;- 7 ln(x))) 7 _\ , K J> = 1- ( 55 ) 

^ ' n=l 

Finally the variables ^ are generated as described above. In this algorithm all events are 
generated with weight 1, provided x is generated exactly according to p(x), for instance 
using the general-purpose tool FOAM. The algorithm is very efficient and fast. 

Numerical results from the corresponding MC program are shown in fig. |3J Plot (A) 
in this figure shows the distribution which is the integrand of eq. (J5U)) from the MC run. 
The analytical result superimposed on the same plot is indistinguishable from the MC 
result. In the next plot (B) we show the ratio of the two distributions, MC and analytical. 
They agree within a very small statistical error of order ~ 10~ 4 . In the next two plots, 
(a) and (b), we see that the new algorithm reproduces perfectly well the semi-exclusive 
distributions of the same two plots in fig. 01 The multiplicity distribution in the plot (c) 
is also well reproduced. Plot (d) shows the MC weight distribution. 

3.3.3 Prototype benchmark type II. a, pure bremsstrahlung 

The two toy MCs from previous sections should be regarded as introductory exercises (and 
numerical benchmarks) for the next step, in which we shall elaborate on the constrained 
Markovian MC solution with x-tagging of the type II. a. We shall restrict ourselves to 
pure bremsstrahlung from the gluon or quark line, without using the multibranching to 
isolate Pq G ~ 1/z. The corresponding MC prototype we name BremsP. The purpose of 
that is threefold: (a) to measure the MC efficiency of this class of the MC algorithms, (b) 
to provide a cross-check for the more sophisticated prototype MC with multibranching 
for the bremsstrahlung from the gluon line, which will be developed in the next section, 
(c) to compare it with the other constrained Markovian MC prototypes for the pure 
bremsstrahlung. 

The starting point for the construction of the algorithm is eq. (jl4J) . Its simplified 

u It is done by using the simple/universal method of inverting cumulative distribution. 
12 Lct us note that a similar Bessel-type distribution of the number of emissions is used by Kharraziha 
and Lonnblad in the event generator based on the Linked Dipole Chain model ^S] ■ 
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Figure 5: Results from BremsP for the gluon emitter, k = G. 

version, restricted to the pure bremsstrahlung case, is the following: 

T 1 

/oo n „ „ 

dx H k (x) A k x~ 1+ v e~( T -^ Rk >< dTi dz * ? kk(zi) 

n=0 i=l J „ 



n-i x 



where 



W?(x,x ) 



xW W k D (x,x o ) J 
D k (r ,x ) f x 1 ~ n 



A--1+V 1 = < ■ 
■rt-k.Xn 



w e =Q\ n 



\. k x lx 

Neglecting W e W D we can perform a z-integration and ra-summations 

i 



Zj — x | . 



dx H k {x) A k x~ 1+r] e {T - To)( - Rkk+n ^ x)) , Q kk (x) = dz ?® (z). (58) 



o 



The emission multiplicity distribution is Poissonian: 

P n {x) = i-e- A ^A(x) ra , \{x) = (r - r )Q kk (z) 



(56) 



(57) 



(59) 
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Figure 6: Results from BremsP for quark emitter, k = Q. 



and we may generate it together with the T{ variables, much as in the Markovian case, ex- 
cept that the average multiplicity (forward leap in Markovian random walk) now depends 
on x (in the unconstrained Markovian it was constant). The variable x is generated as a 
first variable using FOAM then n and finally Zi G (x, 1) exactly according to Tf^z). 

A few comments on the form factor are in order here. The part (r — r )fifefc(x) is clearly 
coming from the real emission and, for instance, will be different if we generate according 
to an approximate 7f k (z)] see later in this section. The part —(r — T )Rkk = —&kk(T, To) is 
a genuine virtual part of the form factor, independent of any details of the MC generation, 
cf. eqs. (P|)-(fTD|). With the usual expansion 



Pik(r,z) =5(1- z)S ik A 



Ti ^—SikBkk H — C%k + D ik (z), 

(l-z)+ z 



(60) 



we obtain 



Rkk = (t - t ) l ^ kk {r, T ) 



a s {t j 



and the real emission form factor is 



7T 



Bkk m A k k 

e 



(61) 



tt k k(x) 



a s {tj 



71 



1 — X 1 

B k k In h Ckk In - + / dz D k kiz) 

e x 



(62) 
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Figure 7: Results from the non-Markovian BremsP (II. a) compared with results of the Marko- 
vian EvolMC. Evolution from 1 GeV to 1 TeV due to multiple gluon emission from the gluon 
emitter line (upper plots) and quark emitter line (lower plots). Starting D k (r ,z) as in the 
realistic proton. Nf = 3 in the virtual form factor. 



where 

-l 



dzD GG (z) =2C A (^- l ^ + x(2- l -x + ^x 2 yj (63) 
and ^ 

dz D qq (z) = C F (^-^ + x + ±x 2 ) . (64) 

In fig. we show type II. a MC results for the same semi-exclusive distributions as 
previously, using realistic gluon distribution Dg(t ,x) = Cx~ 8 (l — x) 5 , for the gluon- 
strahlung out of the gluon emitter line. As we see, the efficiency of the MC is extremely 
low - the acceptance rate is merely 1.5 x 10~ 5 (note that the weight-0 events are not 
included in fig. Eli). Nevertheless, these results will still be useful to cross-check the more 
efficient algorithm type II. b in the next section. We have investigated what the sources 
of the inefficiency are. As in the previous toy model, the main reason for low efficiency 
is that there are many zero- weight events due to W e . The factor (1 — x) 5 in the gluon 
distribution causes a loss of efficiency of a factor 3. The factor x -0 ' 8 accounts for a mere 
factor 2 in the efficiency loss. It is therefore not urgent to eliminate this efficiency by 
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means of incorporating z~ v factor into Dqg{z)- This possibility we have considered in the 
general discussion on method II. On the other hand, a factor 3 loss in the MC efficiency 
in method II, which is due to the presence of (1 — z) 5 in the gluon SF, looks at first 
sight irreducible. Nonetheless, one may consider modelling this factor using the internal 
rejection loop, because the (1 — z) 5 factor, upon expanding, is a sum of monomials z v and 
the overall normalization can be calculated (with non-MC methods) as a sum over these 
terms. It is not excluded that with some extra effort, the overall efficiency of the method 
II. a could be improved to the level of 10~ 4 . 

Let us now repeat the same exercise for the bremsstrahlung emitted from the quark 
line. In fig.|)]we show the corresponding results (k = q) and the starting quark distribution 
being D q (r , x) = D sea (z) + Du(z) + Dd(z), that is sea plus both valence quarks, taking 
a typical parametrization of the proton parton distribution function at Qo — 1 GeV. 
Strikingly, the overall efficiency is very good; the rejection rate ~ (w)/(w max ) is only 
about 30%! Obviously, without 1/z component in the kernel, the basic algorithm type II is 
quite efficient. It should be remembered that in the actual run P qq (z) is generated exactly 
(i.e. with the help of the internal rejection loop). The fact that the weight distribution 
extends above 1, up to 2.5, is related to the valence component. However, the entire 
weight distribution looks very well for the optional rejection method. 

The overall normalization of this MC is cross-checked with the help of the Markovian 
MC EvolMC of ref. |S]. In the top part of fig. Owe show result of the evolution from 1 GeV 
to 1 TeV in which we restrict ourselves to gluon emission out of the gluon line, taking the 
starting gluon distribution as in the proton. The non-Markovian type II. a MC BremsP 
reproduces the results of the Markovian MC EvolMC within a statistical error of a few per 
cent. The apparent discrepancy at high x values is most likely due to some technical bias 
related to extremely high MC event rejection rate 13 . 

In the low part of fig. [7|we present the analogous comparison of BremsP and EvolMC 
for multiple gluon emission from the quark line. Again, the agreement is quite reasonable, 
this time within a smaller statistical error of ~ 1%. 

As an additional cross-check we also implemented another variant of the II. a type 
constrained MC algorithm BremsP, with the approximate kernels 7 and correcting weight 
applied at the very end of the MC generation. In this case we define 

r 1 

/OO Tl/ n n 

di^(x)V +, e^*'yT] (in dziTf^Zi) , 

n=0 i=l J J ( 65 ) 

n-i x 

xW e W°(x,x )W*, 
where the additional weight is 

13 We did not try to investigate its precise source, because the practical importance of BremsP is limited 
to a test of semi-exclusive distributions, not normalization. 
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Neglecting weights we have 



i 



a,. i dxH k {x) Ak x- 1+ ^ e^- T ^- R ^ +n ^\ Q kk (x) = / dz?g(z), (07) 



where the simplified kernel is defined as 



Pkkit, z) — 5(1 — z)A kk + — r-B kk H — C kk , 

1 — Z 4- 



1 

-( 

z 



(6* 



leading to the following real emission form factor 



Q kk [x) 



as(tA) 



7T 



1 — X 1 

B kk In h C kk In — 

e x 



(69) 



We have checked that the above MC algorithm gives the same quark and gluon distribu- 
tions, as expected. It is also quite interesting to check how strongly the efficiency of the 
MC deteriorates when the additional weight W v is introduced. In the quark case, the 
acceptance rate drops from 0.7 to 0.25, which is not much, while for gluons it drops by a 
factor ~ 10, well below 10 -5 . 

In the next step we will clone the MC subgenerator of type II, which generates brems- 
strahlung from the quark line according to simplified P® = 2CV/ (1— z) and from the gluon 
line according to "truncated" simplified Pqq = 2(7^/(1 — z). After that, having tested 
the components at hand, we shall introduce the integration over Z using FOAM and for the 
bremsstrahlung from the gluon line we shall combine the Bessel's MC with Pqq = 2Ca/z 
with the above MC for = 2C A /(1 - z) and compare resulting distributions with the 
Markovian benchmark of fig. [7| This will close the most important first step in making a 
prototype MC according to method II. b. 



3.3.4 Constrained MC type II. b for pure bremsstrahlung 

In the following we implement the algorithm II. b in the case of pure bremsstrahlung from 
the gluon or quark line. In this particular case, the master formula of eq. (J39|) for the early 
stage MC (obtained from eq. (JHEj) by neglecting the MC weight) has only one variable 
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and z^ in = x/Z^. It takes the following simplified form: 

1 

o k = fdxH^x) A k x~ 1+V J dZ<» e~^o)K df (r, z^M J^rff ( T ,Z«|r ), 

X 

co n' T „ n' \ 

ft4w=e- ( ^E n / d <n n / ^^(o 



n'=0 m=l , m=l 



eff 
min 



exp {-{t - roX, + (r - r ) jf tfetff^*) j 

2 min 



<(r,Z«|r ) = 5(l-Z( 1 )) (7Q) 
dg(r,ZW|r ) = e-^ R oo!^6(Z^ - 1) + 



n« 1 



IV n IV n IV 

+ E II / d^ll dz' m n^'js(z^-l[z 

n'=l m=l , m=l n m=0 

T m-1 U 

= ^'( (r - ro) ^' Z(1) ) 

= e - {T -^ R oo{5(Z^ - l) + d R (r,Z^\r )), 

dg(r, ZW|r ) = ^y(r - r )^ G Fi(2; -(r - r )i^ G ln(Z«)) . 
The integral proportional to 5(1 — Z^) has to be treated separately 14 : 

a q = f dx H q {x) A q x- 1+r > e -(— o)^ df(r,x\r ), 
Jo 

i 

*a = f Q dx H G (x) A G x-^ j ^e'^^c (r, ^\r ) e^^r, Z^\r ) 

X 

1 

+ I ix H G(X ) ^ '°4\r,z Me -^° U. 



(71) 

The distribution of the variables x and Z = Z^ for the general-purpose simulator FOAM 

14 The need of treating the <5-part separately will be more annoying in the general case, with several 
gluon emitter bremsstrahlung segments, because this causes proliferation of the separate MC branches 
with different distributions, adding a lot of code, difficult to write and debug. 
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are given by the integrands in the integrals: 
1 1 

a q = dx dZ H q (x) A q x~ l+V exp 



x 

i i 

dx I ^ H G (x) A G x- 1+ " 



r-r )(R'+Rf- I dz?Z A (z 







exp 



(r - r ) ( R' G + R GG + R» G -J dz9 GG (z] 

x/Z 



d G (r,Z\r ) 



r 

+ J dx H G (x) A G x- l+ ^ exp -(t-t )(r g + R% g + R£ g - j dz? GG (z)^j 



(72) 



Keeping in mind that Ri 3 = 0, we recover in eq. (J72|) the complete virtual form factor Rk 



R' k + Rt + R% k = R k = (r- r ) ^ k (r, r ) 
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Bkk In A kk 

e 



(73) 



see eq. (16 Finally we arrive at the following expression: 



i i 



o q = dx dZ H q (x) A qX - 1+TI e 



x 

a (a) a (6) = 

1 1 



a G = jdxj^ H (x) A G x~^ e-^( R -- n o^)d G (r,Z\r ) 



(74) 



x 
1 

+ dx H t 



G 



(x) A G X~ l+rt e -(T-T Q ){R. G -%{x)) 



— °G + °G ■ 

The MC algorithm of type II. b for generating single (weighted) MC event consists of 
the following steps: 

— (X) 

1. Generate a branch index X = a,b according to a probability proportional to a k ; 
FOAM does that efficiently 

2. For given X generate variables x and Z or only x according to the integrand of the 

PO. 



corresponding integral a k , also done by FOAM 
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Backward evolution 



Forward evolution 
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Figure 8: Scheme of the II. b MC algorithm in x space before relabelling, for the case of pure 
bremsstrahlung. x = x /(ZTli z i) faults from Z and Yl z i- The zero-weight case of xq > 1 
is also indicated. 

3. In the case X = a generate two emission multiplicities n' and n*, the first one 
according to the Poisson distribution with (n') = (r — t )Q^(x/Z) and the other 
one according to the Bessel-type distribution with A = (r — t q )Rq G ln(l/Z) (as in 
the toy models). 

4. Knowing the multiplicities, generate the variables (r/, z^), i — 1, . . . , n' and (r*, z*), j = 
1, . . . , n*, using methods described earlier. 

5. Relabel the emission vertices, guided by the order of the r variables. 

6. Calculate the final MC weight, the same as was neglected at the early stage of 
generating "phase-space" variables. 

The above algorithm is also illustrated schematically in fig. El in the x-space, before the 
relabelling. Arrows help to understand the order of the reconstruction of all x variables 
out of z variables. 

In fig. El we show numerical results for the II. b prototype MC for the same semi- 
exclusive x- and r-distributions as previously. MC results coincide very well with these 
from BremsP in figs. El This is a highly non-trivial result, having in mind sophistication 
of the algorithm II. b. Let us stress that the above agreement cannot be obtained without 
a correct relabelling procedure being performed in the final stage of the algorithm II. b 15 . 

The generation time of an event (before any rejections) is similar for both algorithms 
II. a and II. b. Therefore the acceptance, i.e. the ratio of the average to maximum weight, 
is a good measure of the overall efficiency of the algorithms. The acceptance for the new 
algorithm type II. b, as read from the weight distribution in fig. El is 9.6 x 10~ 3 . This is a 
little bit worse than expected; it is, however, fully satisfactory - it is better by a factor of 
500 than the efficiency 2 x 10~ 5 for the solution II. a (without multibranching), see fig. El 

15 We have checked this fact numerically in a separate MC exercise. 
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Figure 9: Results form Genllb MC prototype for gluon emitter, k = G. 



Using algorithm II. a as a guide, one may argue that the 9.6 x 10~ 3 efficiency can still 
be improved by a factor of 2 by including exactly the Xq °' 8 factor. Another factor of 3 
could be obtained by performing a modelling of the (1 — Xq) 5 distribution in an internal 
rejection loop of the algorithm. In this way the overall efficiency may go up to the level 
of 5%. 

Plots in fig. El show tests of exclusive distributions and efficiency, but not the overall 
normalization. A strong test of the overall normalization of the algorithm II. b is shown 
in fig. where high statistics (~ 4 x 10 9 events) results of the II. b MC are compared 
with those of the forward Markovian MC EvolMC of ref. [8j 16 . The agreement is reached 
within a statistical error of about 0.1% for x < 0.01 and of 0.3% for x < 0.1. For higher 
x, in spite of the extreme smallness of Dg{x) (over 9 orders of magnitude), the agreement 
holds perfectly well within the statistical errors. 

At this point we may state that our method II. b of solving the constrained MC really 
works in practice and is reasonably efficient. 

We want to stress that it would not be possible to reach this conclusion without 
constructing and testing the explicit prototype of the algorithm type II. a, and other 

16 Results of EvolMC were in turn cross-checked very precisely with the results of two non-MC evolution 
programs; see ref. [Hj- 
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Figure 10: Comparison of II. b MC Genllb with the Markovian MC EvolMC for pure gluon- 
strahlung k = G, Dg(t , z) as in proton. Nf = 0. 



auxiliary MC exercises, as we did in this work. 

In the above numerical exercises we have restricted ourselves to the LL case, with 
Nf = 3 massless quarks. The QCD evolution kernels are unique and well known, and we 
therefore skip their explicit definition. The running constant ats(t) = 2ir/(f3 (t — hiA )) 
was used with A = 0.245748338. The following starting values of the parton distributions 
in proton at Qq = 1 GeV were used in all our numerical exercises: 



xD G (Q ,x) 
xD q (Qo,x) 
xD q (Q ,x) 

^-Dsea(Qcb x ) 

xD 2u (Qo,x) 
xD d (Q ,x) 



1.9083594473 • aT a2 (l - x) 5 ' , 
0.5 ■ xD sea (x) + xD 2u {x), 
0.5 ■ xD sea (x) + xD d (x), 
0.6733449216 • x" a2 (l - x) 7 , 



(75) 



2.1875000000 ■ x 
1.2304687500 ■ x 



0.5/ 



0.5/ 



X 

x) 



,3.0 



4.0 
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4 Summary 



In this paper we presented a MC algorithm, which belongs to a new class of MC algorithms 
capable of generating a constrained Markovian evolution of parton distributions according 
to DGLAP evolution equations. Practical numerical implementation is for the moment 
restricted to the pure bremsstrahlung case. Since the algebraic framework is defined for 
the full DGLAP, it is therefore a matter of more programming to extend it to the general 
case. In the presented numerical tests (pure bremsstrahlung) the algorithm has been 
checked against the dedicated forward evolution (Markovian) MC program that we have 
written. We found an agreement at the level of 0.1%. The measured efficiency of the 
constrained MC is found to be quite satisfactory. This work opens the way to a new class 
of MC algorithms with possible applications in the initial-state QCD parton shower MC. 
Furthermore, the Bessel-type distribution of the number of emissions, which forms the 
core of our algorithms, is similar to this obtained from the CCFM approach as shown 
in |18j . This suggests another possible area of applications of our algorithms. 
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APPENDIX 



The technique of the kernel split (multibranching) 

In section 13.2.21 we have shown how to reorganize integration variables in the evolution 
iterative solution, such that in the Monte Carlo integration/simulation algorithm it is 
possible to generate first the chain of flavour indexes (gluon or quark type) and the cor- 
responding evolution time variables r, (i.e. those of the emissions which change flavour), 
and later the other variables corresponding to gluon emissions (no flavour change). 






j_ _!_ 

l-z z 

Figure 11: Graphical representation of the split of the (approximate) gluon emission kernel 
into two parts. 

In the following we shall describe the application of the MC technique of multibranch- 
ing to our problem. In the multibranching one splits the integrand into many positive 
components, chooses randomly one at a time and generates points according to this par- 
ticular component. In the context of the iterative solution of the evolution equation, it is 
worthwhile to apply this technique to the kernel for the transition of the gluon into gluon: 



P GG (z) = 2C A 



+ + z(l - z) 



(76) 



which has two very different singularities \jz and 1/(1 — z). Therefore it is profitable 
in the Monte Carlo to split P GG (z) = P GG (z) + P GG (z) such that P GG (z) ~ 1/(1 - z) 
and Pq G ~ 1/z (see fig. ITTj) . and to generate them separately, applying additional MC 
methods suited to the individual character of each type of singularity 17 . 

Since we already know from section 13.2.21 how to isolate the pure bremsstrahlung 
subintegrals d' k ; see eq. (fT5|) . let us concentrate on one of them: 



d' k (r,x\r Q ) = e - {T - To)Rkk 5(x - 1) + J^l / dr 3 Gfa - r^) / dz 

«=iW/ o { 

' n \ / n 

K i=\ ) ^ i=l 



(77) 



17 One should also take care of the positivity of the two components. For simplicity we would like to 
have Pq G (z) = 2Ca/z. However, in such a case Pa G (z) — Pgg(z) — Pqg( z ) ^ s n °t positive. A possible 
solution is to first simplify Pgg{z) — > Pgg{z) = 2Ca[1/z + 1/(1 — z)], compensating the simplification 
with the MC weight at a later stage, and then to split Pgg{ z ) into two positive components without any 
problem. We shall come back to this point later on. 
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where in reality we are interested in the gluon case Vkk(z) = Vgg{z)- 

Taking advantage of the independence of the kernels on r we can rewrite the above 
equation as follows: 



d' k (r, x\r ) = e -<™»fl» £ iL-pL JJ j dzi 

n=0 ' i=l Q ^g) 

n / n \ 



x I I I /-/...( :,) + /-;./( <)[ .c - I I ,. 

i=l x i=i 

o 

where more compact notation is achieved by defining fj = 1. Note that at this stage we 

i=i 

made certain important short-cuts, because we have integrated over r. This simplifies the 
argument but makes it questionable in view of certain important claims concerning the 
final distribution in the space of (n; Ti, zi, t%, Z2, ■ T n , z n ), which we are going to make at 
the end of this appendix. We shall therefore refine our proof later on, showing how to 
proceed for the distributions with unintegrated r's. 

□ o 



Figure 12: Reorganization of the multiple gluon emission leading to multibranching. 

Let us now reorganize the overall sum as follows (see also schematic illustration in 

fig. inn, 

i i 

oo oo / \ rai _|_ n2 ni „ ri2 „ 

4(t,.m = e -<->«» £ £ (T "t, nt*»nt ^ 

ni=0n 2 =0 i=l Q j=l Q ^ 

ni n.2 / ni ni \ 

i=l j'=l ^ i=l j=l 

where the two sums take care of the two kernel components. We can now factorize the 
whole integral as a convolution of the two integrals, each of them corresponding to one 
component of the kernel: 



4( r ^l r o) = / 

JO 

oo / \n n / n \ 

d' k x (r, x \r ) = e-(-o)^ £ J] 5 [ x ~ U ** ) > X = A > B - 

n=0 ' i=l ^ j=l ' 



dz A I dz B 5{x - z A z B ) d' k A (r, z a \t ) d fB (r, z B \r ), 

(80) 
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The functions R x are constrained only by 

X 

For example in some cases they may be defined as 

-l-s 



J2& = / dzzVtCiz)- (82) 
Jo 

We may restore the ordered evolution time integrals 
However, it should be remembered that the variables Tj and z% are not exactly the same as 



oo n p n / n \ 

e -(r-roM k J- H / dr j Q{r j -T j - 1 )\{vg x {z{)8L-\{zX (83) 

n=0 7=1 ^ i=l ^ i'=l ' 
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Figure 13: Further reorganization of the multiple gluon emission in the multibranching. 



in the original integral (in spite of the same notation) but they are related by means of a 
"relabelling" procedure described later in this appendix. The above algebra is represented 
schematically in fig. ED 

It is now possible to implement the integral of eq. (|50|l as a pair of two independent 
"parallel Markovian processes", both starting at r and stopping at r. The first one would 
have decay constant R^ k and variable Zi generated according to V^(zi), yielding emission 
multiplicity ni at the stopping point, while the second one would have its decay constant 
R^ k , variables Zj generated according to V k 3 k B (zj), and the emission multiplicity n-i- 

It is important to understand that at the very end the two sets of Zi,Ti, i = 1, ...,n\ 
and Zj,Tj, j = 1, ...,n 2 can be merged, forgetting from which parallel generation branch 
they originate. Merging is done simply by creating a common list of ordered variables T; 
and renaming/reordering Zi variables in exactly the same way. Such relabelling procedure 

will undo the procedure of combining together the (^ ni ^ 2 j terms done in eq. (|79|). The 
relabelling procedure is illustrated schematically in fig. The resulting Zi, I = 1, ...,ni + 
n 2 will be then distributed according to the product 

711+T12 

II [Vkk\zi)+VT(zi)} (84) 
i=i 

Moreover, also the total multiplicity n = n\ + n 2 and the evolution times 77, I = 1, ...,n 
will be distributed as if they were coming from the corresponding single Markovian MC. 
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Figure 14: Illustration of "relabelling". The actual generation is done in two steps: First, 
for each of the two branches (squares and circles) the ordered r's are generated separately 
and independently in the entire r-range. Next, {r^Zj) are relabelled according to a common 
ordering in r. Only after such a relabelling is x constructed: x = YYj=i z j- 



Actually, the reader may be concerned that the above claim is not really founded on 
a solid derivation because we have excluded the r space in the binomial decomposition 
after integrating over r's at an early stage of derivation, while we are now making state- 
ments on the distribution in the full space n, Z\, r 2 , z 2 , T n , z n . We need clearly to refine 
our derivation keeping the r-space alive. The full derivation involves non-trivial combi- 
natorics, and here we shall only give a sketch on the necessary reasoning. Consider the 
expression with three kernels 

dT 1 dT 2 dr 3 e 321 J dz 1 dz 2 dz 3 {A{3) + B{3)){A(2) + 5(2))(A(1) + 5(1)), (85) 

where we abbreviate: 6 32 i = 0(r 3 — t 2 )6(t 2 — T\) and A(i) = ^^(^j), B(j) = "Pkk( z j)- ^ 
is decomposed as follows 

dr 1 dT 2 dT 3 Q 3 2i / dzxdz 2 dz 3 [A(3)A(2)A(i) 



+ B(3)A(2)A{1) + A(3)B(2)A(1) + A(3)A(2)B(1) ( 86 ) 
+ B(3)B(2)A(1) + B(3)A(2)B(1) + A(3)B(2)B(1) 
+ 5(3)5(2)5(1)] 

Each of the four groups in four lines is now transformed separately into a single factor with 
different ordering pattern of the r variables. For instance the second line we transform 
explicitly as follows: 



dndr^Q^ J dz 1 dz 2 dz 3 [B(3)A(2)A(l) + A{3)B{2)A{1) + A{3)A{2)B{1)\ 
dridT 2 dTrdz 1 dz 2 dzvQv 21 B(l')A(2)A(l) 



(87) 

+ / dTidT V dT 3 dzidz v dz 3 <c) 3V iA(3)B(l r )A{l) 



+ / dr v dr 2 dr 3 dz v dz 2 dz 3 e 32V A(3)A(2)B(l')} 
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where we essentially renamed both z's and r's sitting in the I?-factor. The same can be 
done for variables in the A-factors: 



dT'dT'dT' l dz 9 1 dz'dz' 1 Q l/2 . 1 .B{l')A{2')A{r) 
+ J dTldT[dT*dzldz 1 dzlQ 2 . vl .A{2 9 )B{l')A{l m ) 
+ J dT[dT^dr'dz[dz , 1 dz , 2 e 2 . 1 . 1/ A(2')A(r)B(l')]. 

Let us summarize explicitly the relabelling of the variables that has been done above: 



for t[ > r* > t* 
for r* > t[ > t* 
for r* > r" > r[ 



T\ = t[, t 2 = t', n = r*, zi = z[, z 2 = z 2 , zi = zl, 
Ti = r', r 2 = t[, Ti = t', Zi — z*, z 2 = z[, Zi — zl, 
Ti = r 2 , r 2 = r', n = t[, zi = z 2 , z 2 = zl, Zi = z[, 



(89) 



Now, we may pull out the kernels and combine the 0-functions 

J dT'dT 2 dridzldz 2 dz[B(l , )A(2')A{l')[e V2 . 1 . + Q 2 . V1 . + e 2n . 
dT'dT'dr[ dzldz' 2 dz[ e 2 . 1 .A(2*)A(l*) e[B(l'), 



(90) 



where Or = 1 and 2 »i« = 0(t* — t'). The above tedious relabelling of z's and r's and 
recombining of 0's into product of two independent ones can be done for any number of 
kernels. The net result is an interesting identity: 



[ [ / dr { d Zi (A(zi) + B(zi))w(T,z) 
i=i^n-i Jo 

n / m „i „i \ /n-m pi i>i \ /g^\ 

= E(n < / T . dr' J o dz' A(z')j fj|| d^ jf dz[ A{z'^\ 

wHT'.TVirV'z'.z')), 

where w(t, z) is an arbitrary "test function" ensuring that eq. (j91)l is indeed a differential 
identity, and not an obvious multiplication rule of exponential functions. The mapping 
(relabelling) r» = Tj(r', r') and z% = Zi(r'; t' , z', z') is nothing more than a permutation 
of the integration variables, which is "guided" by the ordering of the r variables, much as 
in the explicit example above. Note that the above identity is still valid if the integrand 
involves any additional factor, symmetric with respect to the permutation of the integrand 
variables Ti and z^. 

The above formula is a kind of generalization of the Newton binomial formula in 
which an n-dimension simplex in r variables is decomposed into a sum over the Cartesian 
product of the two simplexes in ni and n — % dimensions. From this exercise it is also 
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clear that this identity implicitly involves a relabelling of z variables depending on the 
ordering of r variables. This is exactly what we have to do in the Monte Carlo if we 
generate A(z) and B(z) independently, but we want to have the distribution A(z) + B(z) 
at the end of the algorithm. Note that a similar MC procedure with relabelling of the 
integration variables was done in the context of the ISR and FSR photon radiation in the 
YFS3 algorithm of KKMC, before adding the ISR-FSR interference [H]. 
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